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Abstract 

The dynamics of Boolean networks (BN) with quenched disorder and thermal noise is studied via 
the generating functional method. A general formulation, suitable for BN with any distribution of 
Boolean functions, is developed. It provides exact solutions and insight into the evolution of order 
parameters and properties of the stationary states, which are inaccessible via existing methodology. 
We identify cases where the commonly used annealed approximation is valid and others where it 
breaks down. Broader links between BN and general Boolean formulas are highlighted. 
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In his seminal work [l| Kauffman introduced a very simple dynamical model of biological 
gene-regulatory networks. The state of each gene was modeled by an ON/OFF variable, 
interacting with other genes via a coupling Boolean function which determines the state of 
a gene at the next time-step. There are N such genes (sites) in the network and each gene 
is influenced by exactly k other genes from the same network. In Kauffman's approach, the 
networks are constructed in a random manner by choosing Boolean functions from the set 
of all 2^* functions of k inputs and by connecting the inputs of each function to the genes 
randomly selected from the set 1, .., A^; Boolean functions and connections are fixed for all 
subsequent time-steps {quenched variables). The evolution of a such dynamical system is 
deterministic and since the number of states is finite (2^) the system is driven to a periodic- 
orbit attractor. 

It was argued |l| that, despite its simplicity this model, also known as Random Boolean 
network (RBN) or Kauffman net, is of relevance to the understanding of biological systems 
a„d ^ been .ud.d p„.„i,y fo. .hi. ..on « . RBN belong to a la... Ca. of Boolean 
networks, the N-k model of A^-variable dynamical systems with a discrete state-space and 
/c- variable interactions, that exhibits a rich dynamical behavior ^, Q]. The N-k model is very 
versatile and has found its use in the modeling of genetic networks [sl, neural networks jo], 
social networks [\ and in many other branches of science js, 4|. 

For over two decades the annealed approximation [s^ has proved to be a valuable tool in 
the analysis of large scale Boolean networks (A^ — ?■ oo) as it allows one to predict the time 
evolution of network activity (proportion of ON/OFF states) and Hamming distance (the 
difference between the states of two networks of identical topology) order parameters. The 
latter was used to predict a phase transition at A; = 2 in RBN. The main assumption in 
this method is to ignore the fact that both Boolean functions types and random connec- 
tions in a Boolean network are quenched variables and enables one to resample them at each 
time-step. This allows one to ignore the correlations among input-variables, which simplifies 
an analytical treatment significantly. It was .hownBH tha; the annealed app^o^iniallon 
indeed gives a correct result for the Hamming distance order parameter in RBN, but the 
broad validity of the annealed approximation to general networks of this type has remained 
an open problem {ll]. Remarkably, the annealed approximation provides accurate activity 
and Hamming distance results for many other Boolean models with quenched disorder but 
cannot compute correlation functions, used in studying memory effects, due to the repeated 
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resampling at different time steps that makes the various quenched systems indistinguish- 
able. Furthermore, there are models 1^ that have very strong memory effects in specific 
regimes, where the annealed approximation is no longer valid. 

In this Letter, we study the dynamics of the N-k model with quenched disorder and 
thermal noise using the generating functional analysis (GFA), an established method for 
studying physical systems of this type the analysis is general and covers a large class 
of recurrent Boolean networks and related models. We show that results for the Hamming 
distance and network activity obtained via the quenched and annealed approaches, for the 
N-k model, are identical. In addition, stationary solutions of Hamming distance and two- 
time autocorrelation function (inaccessible via the annealed approximation) coincide, giving 
insight into the uniform mapping of states within the basin of attraction onto the stationary 
states. In the presence of noise, we show that above some noise level the system is always 
ergodic and explore the possibility of spin-glass phase jl4| below this level. Finally, we show 
that our theory can be used to study the dynamics of models with strong memory effects. 

The model considered is an A^- variable recurrent Boolean network with the parallel update 
rule 

S,{t+l) = ai{S,,{t),...,Si,{t)), (1) 
where Si(t) G {—1, 1} and : {—1, l}'^ —J- {—1, 1} is a Boolean function of exactly k inputs. 



We assume that the thermal noise can flip the output of a function with probability p [15|. 
The function at site i and time-step t + 1 operates in a stochastic manner according to the 
microscopic law 

P^XS^{t+l)\S.At),..,SM) (2) 

g/35,{t+l)a,(5,^{t),..,S,Jt)) 



2 cosh (3ai{Si^{t), .., 5'i^(t)) 

where the inverse temperature /3 = 1/T relates to the noise parameter p via tanh /3 = 1 — 
2p. The function-output Si{t + 1) is completely random/deterministic when /3 — )■ O/oo, 
respectively. Given the state of the network S{t) G {—1, 1}^ at time t the functions at time 
t+1 are independent of each other. This suggests that the probability of the microscopic path 
S{0) — 7- ■ ■ ■ — S{tmax) IS a product of ([2j) over sites and time steps. The joint probability 
of microscopic states in two systems of identical topology but subject to different thermal 
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noise is 



P[{5(t)};{5(t)}] = P(5(0),5(0)) (3) 

imax 1 

X Yl PiS{t+l)\S{t))P{S{t + l)\S{t)) where, 



t=o 



p{S{t+l)\s{t))=u^,PaAs^{t+l)\s,At),■■,s.M■ 

The quenched disorder in our model arises from the random samphng of connections and 
Boolean functions generated by selecting the i-th function and sampling exactly k indices, 
{ii, .., ik}, uniformly from the set of all possible indices. Boolean functions {ai} are sampled 
randomly and independently from the set G of A;-ary Boolean functions. To analyze the 
typical properties of the system via the generating functional method one defines 



where (. . .) denotes the average generated by ([3]). The generating function (jlj) is used 
to compute moments of ([3]) by taking partial derivatives with respect to the generating 
fields {ilj,{t),iijj{s)}, e.g. {Si{t)Sj{s)) = - lim^^^ ^j^^^-— FfV?; -0]. We assume that 

the system becomes self- averaging for -^^ oo 13] and compute r[i/>;'0], where is the 



disorder average; this gives rise to the macroscopic observables 

N , N 



i=l i=l 



1 ^ 



N- 

1=1 



where m{t) is the network activity (or magnetization 16|), C{t,s) is the correlation between 
two states of the same network and Cuit) (related to the Hamming distance d{t) via d{t) = 
|(1 — Ci2(t))) is the overlap between two copies of the same network. 

Averaging (jl]) over the disorder leads to the saddle-point integral r[. . .] = 

/{dPdP}e^*[^'^] where 



^ = i S)P{S, 5)+log^P(5, 5)e-^^('^"^) . (6) 

5,5 5,5 
For A^— T-oo the averaged generating functional is dominated by the extremum of Func- 
tional variation with respect to the order parameters P(5, 5) provides the saddle-point 
equation 



X { J]P«(5(t+l)|5i(t),..,5fc(t))P«(^(t+l)|5i(t),..,^fc(t)) 



(7) 



The physical meaning of ([7]) relates to the average joint probability of single-spin trajec- 
tories S and S in the two systems P{S, S) = lim^v-s-oo jj J2iLi i^l'^'^ ^^i] ^[^'^ while the 
conjugate order parameter P{S, S) is a constant. Equation ([7]) can be used to compute the 
macroscopic observables (jSD, which evolve in time as follows below, denoting S = {Si, . . . , 5*^) 
and where the magnetization m(t) is computed by ([8]) 



m 



(t+l) = A(m(t))=tanh(/3) 

s j=i 

C(t+1, s + =F„(m(t), m(s), C{t, s)) 



1 + Sjm{t) 
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= tanh2(/3) YJl 
Ci2(t + 1) = F„(m(t),m(t),Ci2(t)), 



l + Sjm{t) + Sjm{s)+SjSjC{t, s) 



{a{S)a{S))a 



(8) 

(9) 
(10) 



Results for the order parameters fl8|)- f[T0|) . in combination with ([7]), suggest that the evolu- 
tion of all many-time single-site correlation functions is driven by the magnetization m{t). 
A similar scenario was observed in recurrent asymmetric neural networks jl^, defined on 
similar topology due to similarity in the equations for m{t) and C{t,s). This is not sur- 
prising since asymmetric neural network is a special case of the N-k model when only linear 
threshold Boolean functions are used. Furthermore, for the stationary solution m = fa{m) 
{m = \imt^oo^(t)) the solutions of q = Fa{m,m,q) (here q = lim;:^fy) lim^^op C(t+T, r) is the 



Edwards- Anderson order parameter, used in disordered systems [IJ] to detect the spin glass 
phase where m = and g 7^ 0) and C12 = Fa{m,m,Ci2) are identical. This suggests that 
there is only one average distance |(1 — g) on the attractor and that all points in the 
basin of attraction uniformly cover the stationary states. 

The annealed model, where connectivities and Boolean functions change at each time 
step ([1]) provides identical results for m and C12 to those of (|8]) and ffTOjl However, 
the annealed correlation function C{t, s) = m{t)m{s), where t>s, is the solution of only 
when networks are constructed from a single function type. 



The annealed result [8| for RBN can be easily recovered from equations (151) - fl 101) using the 
property {a{S))a = for all SgI — 1, l}'^ and {a{S)a{S))a = ,WSj^S where the a average is 
taken over all Boolean functions with equal weight. In this case, the magnetization m{t) =0 
for all t>0 and q = tanh^(/3)(^)'^, corresponding to the stationary solution of (Q, has one 
stable solution q^O for all finite /3 and k. For /3 — > oo (no noise), a transition is observed 
from one stable solution q = 1 for A; < 2 to two solutions q = l (unstable) and q^O (stable) 
for k>2 [sl. 

The unordered paramagnetic phase m = is a fixed point of (|H]) only when 
J2si^('^))a ~ 0- This is a stable and unique solution of ([8]) when tanh/3 < 
|2'=-\/fc((^^-J/2);2'=-V(A;-l)((^':2V2)} = ^(^) ^ ^nd even respectively. To prove 
this [12] we first find a Boolean function x such that /x(^) > fa{^) when m G [0, 1) 
and fa{^) > fxi^) when m G (—1,0]; any function from the set x{S) = sgn[ V^ . Sj] + 
S[0;T.''j=iSjh{S), where 7(^) G {-1,1} and such that E5 '^[0; Ei=i 'S'j]7(^) = [23 satis- 
fies these properties. Secondly, we show that m > /x(^) when m G (0, 1) and fxi^n) > ^ 
when m G (—1,0) (/x(0) = 0) for tanh/3 < 6(fc). Thus, the ordered (ferromagnetic) phase 
m 7^0 is a fixed point of ([8]) (if at all) only for values of f3 and k which satisfy tanh/3 > 6(/c). 
Similar results, for odd k only, have been conjectured using the annealed approximation and 
multiplexing techniques [20]. 

For limt^oo TTT'it) = m, q = is a. fixed point of iff {{J2s '^(^)}^) a ~ ^ which occurs 
only for balanced Boolean functions, with an equal number of il in the output. By similar 
argument to the one used in the previous paragraph we show l3] that for m = the point 
g = is a unique stable solution of ([9]) when tanh^ /3 < b{k). The a-averages in equations 
(|8])-(|9]) can be computed for a uniform distribution over all balanced Boolean functions to 
obtain m{t)=0 for all t>0, which implies g = tanh^(/3) ((^)''(l + 2F:r)- a^r)- The latter 
has only one q = trivial solution for any finite (3 and develops a second q = 1 solution 
only for /3 — 00. Thus, the case of m = 0, q ^0 and finite /3 occurs only (if at all) when 
tanh^ l3 > b{k) and for non-uniform distributions over the balanced Boolean functions. 

The upper bound b{k) computed here for k odd is identical to the one computed for noisy 
Boolean formulas [21I. This is since each site i at time t in our model can be associated 
with the output Si{t) of a k-aij Boolean formula^f depth t which computes a function of 
the associated initial states (a subset of {5*4(0)}) [9]. In the presence of noise, a formula of 
considerable depth (large t) loses all input information for tanh /3 < b{k) and odd k 21j. This 
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suggests that the upper bound b{k), for odd k, is more general and is valid for transitions at 
all m values identifying the point where stationary states depend on the initial states and 
ergodicity breaks. For k even such general threshold is not yet known. 

In model ([1]) the state of site i at time t depends on its states at previous times only indi- 
rectly. In the limit of — oo these dependencies become weak and equation ([7]) factorizes; 
this enables one to calculate the observables of interest (151) - fl 101) . However, in a broad family 



of models 
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231] the state of a site i at a time t + 1 depends directly on its state at time 



t. An exemplar model with strong memory effects used to construct a model of cell-cycle 
regulatory network (A^ = ll) of budding yeast [2^] is of the form 

Si{t + 1) =sgn[hi{t)-2h] + Si{t)6[hi{ty, 2h], (11) 



where hJt) = T.j=iCiA'^ + ^hi^)) ^nd 4. G {-1, 1}. Mean-field theory {N 00) was de- 
rived [12] using the annealed approximation in a variant of this model, where the interac- 
tions {^j} were randomly distributed P{^j = ±1) = 1/2. Significant discrepancies between 



the theory and simulation results has been pointed out [12] for integer h values (in this 
case it is possible that 2h = hi{t)), which was attributed to the presence of strong memory 
effects. Refinements of the annealed approximation method improved the results obtained 
only slightly 
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261 ] but break down in most of the parameter space. 
This model (ITTi) can be easily incorporated into our theoretical framework. The result of 
the GFA ([7j) for this process (with thermal noise) can be obtained by replacing the average 
(■ ■ ■ )a by (■ ■ ■ )^ and the probability function Pa{S{t + l)\Si{t), .., Sk(t)) by 

p^{sit+i)\sity,s,it),..,Skit))= (12) 

QPS{t+l){sgn[h(t)-2h}iS(t)5lh(t);2h]} 

2cosh/3{sgn[h{t)-2h] + S{t)6[h{ty,2h]}' 

where h{t) = J2j=i^ji^ + Sj{t)). In the case of /i G M, the probability function fll2p is 
independent of S{t) and equations fl51)- fni]l) have the same structure as model ffTTl) : the 
a-averages {a{S))a and {a{S)a{S))a are replaced by the averages {sga[h{t) — 2h]) ^ and 
(sgn[/i(t) — 2/1] sgn[/i(tV— 2/1])^ respectively. The equation for m(t) recovers the annealed 



approximation result [l2| (using the relation b(t) = {l+m(t))/2). In Fig. [T] (a,b), we plot our 
analytical predictions for the evolution of m{t) and C(t+t^, t^) against the results of Monte 
Carlo (MC) simulation which use flTT]) . The correlation function C{t+tj^,tj^), in the limit of 
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FIG. 1: (Color online). Evolution of the magnetization (m = m{t)) and correlation (C = C{t + 
twii-w)) functions with time t is governed by ([TT]) . Theoretical results (lines) are plotted against 
the results of MC simulations (symbols) with = 10^. Each MC data-point is averaged over 10 
runs. Error bars are smaller than symbol size. Top: Evolution of m (a) and C (b) for /i G M. In 
(b) we plot C for h = 0.5 and A; = 3. Bottom: Evolution of m (c) and C (d) for /i € Z. In (d) we 
plot C for /i = and k = 2. 

t— j-oOjt^— j-oo, approaches the stationary solution of the overlap function (ITUl) as predicted 
(Fig. [1Kb)). 

The situation is very different when /igZ. Then the magnetization m{t) = P{S)S{t) 
where P{S) is a marginal of ([7]) with Pq— J-P^, is no longer closed as in ([8]), but depends on 
2*-i — 1 macroscopic observables (all magnetizations, all multi-time correlations). Thus the 
number of macroscopic observables that determine the value of Tn{t), or any other function 
computed from ([7]), grows exponentially with time. Annealed approximation results 12| for 
this model when h^Z are only exact up to t < 2 time steps (the equation for 6(1) = (l+m(l))/2 
in our approach and in jl2| are identical) and deviate significantly from the exact solution 
at later times (Fig. [TKc)). A typical evolution of the correlation function C(t+t^,t^) in the 
system (fTTj) when /i G Z is shown in Fig. Wid). 

As BNs are instrumental for our understanding of biological and other complex networks, 
and are directly linked to general Boolean formulas there is a need to develop exact tools of 
greater flexibility that cope with complex networks of variable Boolean functions with strong 
memory effects and emerging correlations. This Letter is the first step in this direction. 
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